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Abstract. We introduce a new and efficient numerical method for mul- 
ticriterion optimal control and single criterion optimal control under in- 
tegral constraints. The approach is based on extending the state space to 
include information on a "budget" remaining to satisfy each constraint; 
the augmented Hamilton- Jacobi-Bellman PDE is then solved numeri- 
cally. The efficiency of our approach hinges on the causality in that 
PDE, i.e., the monotonicity of characteristic curves in one of the newly 
added dimensions. A semi-Lagrangian "marching" method is used to 
approximate the discontinuous viscosity solution efficiently. We com- 
pare this to a recently introduced "weighted sum" based algorithm for 
the same problem [25]. We illustrate our method using examples from 
flight path planning and robotic navigation in the presence of friendly 
and adversarial observers. 



Section 1. Introduction. 

In the continuous setting, deterministic optimal control problems are often 
studied from the point of view of dynamic programming; see, e.g., [I], [8]. 
A choice of the particular control a(t) determines the trajectory y{t) in the 
space of system-states C R n . A running cost K is integrated along that 
trajectory and the terminal cost q is added, yielding the total cost associated 
with this control. A value function u, describing the minimum cost to pay 
starting from each system state, is shown to be the unique viscosity solution 
of the corresponding Hamilton-Jacobi PDE. Once the value function has 
been computed, it can be used to approximate optimal feedback control. 
We provide an overview of this classic approach in section [2j 

However, in realistic applications practitioners usually need to optimize by 
many different criteria simultaneously. For example, given a vehicle start- 
ing at a; £ 11 and trying to "optimally" reach some target T, the above 
framework allows to find the most fuel efficient trajectories and the fastest 
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trajectories, but these generally will not be the same. A natural first step 
is to compute the total time taken along the most fuel-efficient trajectory 
and the total amount of fuel needed to follow the fastest trajectory. Com- 
putational efficiency requires a method for computing this simultaneously 
for all starting states x. A PDE-based approach for this task is described 
in section 13.11 

This, however, does not yield answers to two more practical questions: 
what is the fastest trajectory from x to T without using more than the 
specified amount of fuel? Alternatively, what is the most fuel-efficient tra- 
jectory from x, provided the vehicle has to reach T by the specified time? 
We will refer to such trajectories as constrained- optimal. 

One approach to this more difficult problem is the Pareto optimization: 
finding a set of trajectories, which are optimal in the sense that no im- 
provement in fuel-efficiency is possible without spending more time (or vice 
versa). This defines a Pareto front - a curve in a time- fuel plane, where 
each point corresponds to time & fuel needed along some Pareto-optimal 
trajectory. This approach is generally computationally expensive, especially 
if a Pareto front has to be found for each starting state x separately. The 
current state of the art for this problem has been developed by Mitchell and 
Sastry [25] and described in section [3T2l Their method is based on the usual 
weighted sum approach to multiobjective optimization [23]. A new running 
cost K is defined as a weighted average of several competing running costs 
Ki's, and the corresponding Hamilton-Jacobi PDE is then solved to obtain 
one point on the Pareto front. The coefficients in the weighted sum are then 
varied and the process is repeated until a solution satisfying all constraints 
is finally found. Aside from the computational cost, the obvious disadvan- 
tage of this approach is that only a convex part of the Pareto front can 
be obtained by weighted sum methods |13| . which may result in selecting 
suboptimal trajectories. In addition, recovering the entire Pareto front for 
each x G O is excessive and unnecessary when the real goal is to solve the 
problem for a fixed list of constraints (e.g., maximum fuel or maximum time 
available) . 

Our own approach (described in section 13. 3D remedies these problems by 
systematically constructing the exact portion of Pareto front relevant to the 
above constraints for all x E simultaneously. Given £1 C R n and r addi- 
tional integral constraints, we accomplish this by solving a single augmented 
partial differential equation on a (r + n)-dimensional domain. Our method 
has two key advantages. First, it does not rely on any assumptions about 
the convexity of Pareto front. Secondly, the PDE we derive has a special 
structure, allowing for a very efficient marching method. Our approach can 
be viewed as a generalization of the classic equivalency of Bolza and Mayer 
problems [8] . The idea of accommodating integral constraints by extending 
the state space is not new. It was previously used by Isaacs to derive the 
properties of constrained-optimal strategies for differential games |21] . More 
recently, it was also used in infinite-horizon control problems by Soravia [37] 
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and Motta & Rampazzo [26J to prove the uniqueness of the (lower semi- 
continuous) viscosity solution to the augmented PDE. However, the above 
works explored the theoretical issues only and, to the best of our knowledge, 
ours is the first practical numerical method based on this approach. In ad- 
dition, we also show the relationship between optimality under constraints 
and Pareto optimality for feasible trajectories. 

The computational efficiency of our method is deeply related to the gen- 
eral difference in numerical methods for time-dependent and static first- 
order equations. In optimal control problems, time-dependent HJB PDEs 
result from finite-horizon problems or problems with time-dependent dy- 
namics and running cost. Static HJB PDEs usually result from exit-time 
or infinite-horizon problems with time-independent (though perhaps time- 
discounted) dynamics and running cost. In the time-dependent case, efficient 
numerical methods are typically based on time-marching. In the static case, 
a naive approach involves iterative solving of a system of discretized equa- 
tions. Several popular approaches were developed precisely to avoid these 
iterations either by space-marching (e.g., [4C H I29 1 I34"]). or by embedding into 
a higher-dimensional time-dependent problem (via Level Set Methods, e.g., 
|27| ) , or by treating one of the spatial directions as if it were time (resulting 
in a "paraxial" approximation; see, e.g., [28]). For reader's convenience, we 
provide a brief overview of these approaches in sections 12.31 and 12.41 Our 
key observation is that the augmented "static" PDE has explicit causality, 
allowing simple marching (similar to time-marching) in the secondary cost 
direction. 

Our semi-Lagrangian method is described in section 13.41 Since the aug- 
mented PDE is solved on a higher-dimensional domain, any restriction of 
that domain leads to substantial computational savings. In section 13.51 we 
explain how this can be accomplished by solving static PDEs in 0, for each 
individual cost. This additional step also yields improved boundary condi- 
tions for the primary value function in R n+r . 

In section |4] we illustrate our method using examples from robotic naviga- 
tion (finding shortest /quickest paths, while avoiding (or seeking) exposure 
to stationary observers) and a test-problem introduced in [25] (planning a 
flight-path for an airplane to minimize the risk of encountering a storm while 
satisfying constraints on fuel consumption) . Finally, in section [5] we discuss 
the limitations of our approach and list several directions for future research. 

Section 2. Single-criterion dynamic programming. 

Subsection 2.1. Exit-time optimal control. To begin, we consider a 
general deterministic exit-time optimal control problem. This is a classic 
problem and our discussion follows the description in [lj. Suppose f2 C R n 
is an open bounded set of all possible "non-terminal" states of the system, 
while dSl is the set of all terminal states. For every starting state x € f2, 
the goal is to find the cheapest way to leave Q,. 
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Suppose A G R m is a compact set of control values, and the set of ad- 
missible controls A consists of all measurable functions a : R i— > A. The 
dynamics of the system is defined by 

y'(t) = f(y(t),a(t)), 

(1) y(0) = x£tt, 

where y(t) is the system state at the time t, x is the initial system state, 
and / : fl x i h R n is the velocity. The exit time associated with this 
control is 

(2) T x , a = min{t G R+, \y(t) G 50}. 

The problem description also includes the running cost K : £l x A ^ R and 
the terminal (non- negative, possibly infinite) cost q : dQ t— > (-R+ ; o U {+oo}). 
This allows to specify the total cost of using the control a(-) starting from 
x: 

J{x,a{-))= / ' K(y(t),a(t))dt + q(y(T x , a )). 
Jo 

We will make the following assumptions throughout the rest of this paper: 
(AO) Function q is lower semi-continuous and mingn q < +oo. 
(Al) Functions / and K are Lipschitz-continuous. 

(A2) There exist constants fci,&2 such that < k\ < K(x,a) < Jt2 for 

Va; G O, a G A. 
(A3) For every x G fi, the scaled vectogram 

V(x) = {f(x,a)/K(x,a) \ a £ A} 

is a strictly convex set, containing the origin in its interior. 

The key idea of dynamic programming is to introduce the value function 
u(x), describing the minimum cost needed to exit O starting from x: 

(3) u{x) = inf J(x,a(-)). 

a(-)eA 

Bellman's optimality principle [6] shows that, for every sufficiently small 
r > 0, 

(4) u(x) = inf { [ T K(y(t),a(t))dt + u(y(r))\ , 

where y(-) is a trajectory corresponding to a chosen control a(-). If u(a;) is 
smooth, a Taylor expansion of the above formula can be used to formally 
show that u is the solution of a static Hamilton-Jacobi-Bellman PDE: 

min {K(x, a) + Vu{x) ■ f(x, a)} = 0, for 

(5) «(#) = q(x), for a; G dQ. 

Unfortunately, a smooth solution to Eqn. [5] might not exist even for smooth 
f,K, q, and 9f2. Generally, this equation has infinitely many weak Lipschitz- 
continuous solutions, but the unique viscosity solution can be defined using 
additional conditions on smooth test functions fT2l [TT] . It is a classic result 
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that the viscosity solution of this PDE coincides with the value function of 
the above control problem. 

Under the above assumptions the value function u(x) is locally Lipschitz- 
continuous on $7, an optimal control a(-) actually exists for every x £ Q 
(i.e., min can be used instead of inf in formula [3|), and the minimizing con- 
trol value a (in equation [5]) is unique wherever \7u(x) is defined jTj. The 
characteristic curves of this PDE are the optimal trajectories of the control 
problem. The points where is undefined are precisely those, where mul- 
tiple characteristics intersect (or, alternatively, the points for which multiple 
optimal trajectories are available). We will let T be a set of all such points. 
By Rademacher's theorem, T has a measure zero in Cl. 

We note that the assumptions (A2-A3) can be relaxed at the cost of 
additional technical details. For example, if V{x) is not convex, the existence 
of an optimal control is not guaranteed even though the value function can 
still be recovered from the PDE ([5]) and there exist suboptimal controls 
whose cost is arbitrarily close to u(x). On the other hand, if V{x) does not 
contain the origin in its interior, then the reachable set 

1Z = {x € f2 | there exists a control leading from x to d£l in finite time} 

need not be the entire O. In that case, d7Z is a free boundary. We refer 
readers to [T] for further details. 

The above framework is also flexible enough to describe the task of opti- 
mally reaching some compact target TcSJ without leaving 0. To do that, 
we can simply pose the problem on a new domain £l new = 0\T, defining the 
new exit cost q new to be infinite on d£l and finite on dT. The assumptions 
(A0-A3) guarantee the continuity of the value function on fi even in the 
presence of state-constraints; k\ > and the fact that the origin is in the 
interior of V{x) yield both Soner's "inward pointing" condition along the 
boundary of the constraint set (as in [36J) and the local controllability (as 
in [2], for example). 

Remark 2.1. The continuity of u on Cl is a more subtle issue requiring 
additional compatibility conditions on q even if that function is continuous; 
otherwise, the boundary conditions are satisfied by the value function "in 
viscosity sense" only fT]. However, due to our very strong controllability 
assumption (A3), the local Lipschitz-continuity of u in the interior is easy 
to show even if q is discontinuous, as allowed by (AO). Without assumption 
(A3) or its equivalent, discontinuous boundary data typically leads to dis- 
continuities in the value function in the interior as well. Such is the case for 
the augmented PDE defined in section 13.31 

Before continuing with the general case we consider two particularly useful 
subclasses of problems. 

Geometric dynamics. Suppose A = {a E R n \ \a\ < 1} and, for all 
x € Q,a G A\{0}, 

K(x, a) > \a\K(x, a/\a\), f(x,a)=f(x,a/\a\)a, and f(x, 0) = 0. 
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Then it is easy to show that the cost of any trajectory is reduced by travers- 
ing it as quickly as possible; i.e., we can redefine A = S n -i = {a € R n \ 
\a\ = 1} without affecting the value function. In that case, the control value 
a is simply our choice for the direction of motion and / is the speed of 
motion in that direction. The equation ([5|) now becomes 

(6) min {K(x, a) + (Vu(x) ■ a) f(x, a)} = 0. 

A further simplification is possible if the speed and cost are isotropic, i.e., 
f(x,a) = f{x) and K(x,a) = K{x). In this case, the minimizing control 
value is a = — Vu(x)/\Vu(x)\ and © reduces to an Eikonal PDE: 

(7) |V«(a;)|/(aO = K(x). 

The characteristics of the Eikonal equation coincide with the gradient lines 
of its viscosity solution. That special property can be used to construct 
particularly efficient numerical methods for this PDE; see section 12, 41 

Time-optimal control. If K(x,a) = 1 for all x and a, then J~(x,a(-)) = 
Tx,a + q{y{Tx,a))- Interpreting q as an exit time penalty, we can restate 
this as a problem of finding a time- optimal control starting from x. The 
PDE © becomes 

(8) max {— Vtt(x) • f(x, a)} = 1. 

We note that the value function for every exit-time optimal control problem 
satisfying assumptions (AO- A3) can be reduced to this time-optimal control 
case by setting K new = 1 and f new (x, a) = f(x, a)/K(x, a); a proof of this 
for the case of geometric dynamics can be found in |41j . 

A combination of these two special cases is attained when K = 1 and 
/ = 1, leading to a PDE |Vu(a;)| = 1. If the boundary condition is q = 0, 
then u{x) is simply the distance from x to dil. 

Subsection 2.2. Semi-Lagrangian discretizations. Suppose X is a grid 
or a mesh on the domain Cl and for simplicity assume that d£l is well- 
discretized by dX. We will use M = \X\ to denote the total number of 
meshpoints in X. 

A natural first-order accurate semi-Lagrangian discretization of equa- 
tion ([5]) is obtained by assuming that the control value is held fixed for 
some small time r. If U(x) is the approximation of the value function at 
the mesh point x € X, then the optimality principle yields the following 
system: 

U(x) = mm {tK(x, a) + U (x + rf(x,a))} , Va; G X\dQ; 
aeA 

(9) 

U(x) = q(x), Vx G dX. 

Since x a = x + rf(x,a) usually is not a gridpoint, a first-order accurate 
reconstruction is needed to approximate U(x) based on values of U at nearby 
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Discretizations of this type 



were introduced by Falcone for 



time-discounted infinite-horizon problems, which can be related to the above 
exit-time problem using the Kruzhkov transform [171 [IB] . 

In case of geometric dynamics, another natural discretization is based on 
the assumption that the direction of motion is held constant until reaching 
a boundary of a simplex. For notational simplicity, suppose that n = 2 and 
S(x) is the set of all triangles in X with a vertex at x. For every s E S(x), 
denote its other vertices as x s \ and x s2 . Suppose x a = x + r a f(x, a)a lies 
on the segment x s ix s2 . Let E = {£ = (£i, £2) | £1 + £2 = 1 and £1, £2 > 0} . 
Then x a = £ix s % + for some ^GS and U(x a ) = £,iU(x sl ) + £ 2 U(x s2 ). 
Alternatively, given x^ = i\x s \ + ^ 2 x s2 , we can define = {x^ — x)/\x^ — x\ 
and t c = \x^ — x\/f(x, ag). This yields the following system of discretized 
equations: 



U(x) = min min{T^K(x,a^) + (£iU(x sl ) + &U(x s2 ))} , Va; e X\9fi; 



Discretizations of this type were used by Gonzales and Rofman in [18] . Both 
discretizations ([9]) and (|10p were also earlier derived by Kushner and Dupuis 
approximating continuous optimal control by controlled Markov processes 
[23] . In the appendix of [35] we demonstrated that (fTUj) is also equivalent to 
a first-order upwind finite difference approximation on the same mesh X. 

Subsection 2.3. Causality, dimensionality & computational efficiency. 

We note that both of the above discretizations result in a system of M 
non- linear coupled equations. Finding a numerical solution to that system 
efficiently is an important practical problem in itself. 

Suppose all meshpoints in X are ordered x\, . . . , xm and U = (Ui, . . . , Um) 
is a vector of corresponding approximate values. The above discretized equa- 
tions can be formally written as U = J~{U) and one simple approach is to 
recover U by fixed-point iterations, i.e., U k+1 = J-{U k ), where U° is an 
appropriate initial guess for U. This procedure is generally quite inefficient 
since each iteration has a O(M) computational cost and a number of itera- 
tions needed is at least 0{M l / n ). 

A Gauss-Seidel relaxation of the above is a typical practical modification, 
where the entries of U k+l are computed sequentially and the new values are 
used as soon as they become available: U k+1 = Ti(U k+1 , . . . , U^, U k , . . . , Ujfa). 
The number of iterations required to converge will now heavily depend on 
the PDE, the particular discretization and the ordering of the meshpoints. 
We will say that a discretization is causal if there exists an ordering of 



If x is close to dQ, it is possible that x a Cl. This can be handled by either enlarging 
X to cover some neighborhood of 57 (see, e.g., [4]) or by decreasing r in such cases to make 
x a 6 9Q. The latter strategy is also adopted in our implementation of the method in 
section 13.41 



U(x) 



ses(x) £es 
q(x), \ 



Vx e dX. 



(10) 
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meshpoints such that the convergence is attained after a single Gauss-Seidel 
iteration. For example, if the dynamics of the control problem is such that 
one of the state components (say, yi) is monotone increasing along any fea- 
sible trajectory y(t), then ordering all the meshpoints in ascending order 
by x\ would guarantee that only one Gauss-Seidel iteration is needed. Such 
ordering is analogous to a simple time-marching (from the past into the fu- 
ture) used with discretizations of time-dependent PDEs (e.g., in recovering 
the value function for fixed- horizon optimal control problems). If a causal 
ordering is explicitly known a priori (as in the above example), we refer to 
such discretizations as explicitly causal. 

Explicit causality is a desirable property since it results in computational 
efficiency. Suppose u(x) is a viscosity solution of some boundary value 
problem. Historically, two approaches have sought to capitalize on explicit 
causality in solving more general static PDEs. First, it is often possible to 
formulate a different time-dependent PDE on the same domain Q so that 
its viscosity solution <p is related to u as follows: u(x) = t <fr(x,t) = 0. 

The PDE for <j> is then solved by explicit time-marching; moreover, since 
only the zero level set of 4> is relevant for approximating u, the Level Set 
Methods are applicable e.g., see [27], [30]. Aside from increasing the di- 
mensionality of the computational domain, this approach is also subject to 
additional CFL-type stability conditions, which often impact the efficiency 
substantially. Alternatively, in some applications (e.g., seismic imaging) a 
"paraxial" approximation results from assuming that all optimal trajectories 
must be monotone in one of the state components (say, y\) even if the same 
is not true for all feasible trajectories. This leads to a time-like marching 
in x\ direction (essentially solving a time-dependent PDE on an (n — 1)- 
dimensional domain); see, e.g., [28]. This method is certainly computation- 
ally efficient, but if the assumption on monotonicity of optimal trajectories 
is incorrect, it does not converge to the solution of the original PDE. 

Subsection 2.4. Efficient methods for non-explicitly causal prob- 
lems. Finding the right ordering without increasing the dimensionality has 
been the subject of much work in the last fifteen years. This task can be 
challenging for discretizations which are causal, but not explicitly causal. 

For the isotropic control problems and the Eikonal PDE (J7J), an equiv- 
alent of discretization (I10p on a Cartesian grid is monotone- causal in the 
sense that U{ cannot depend on Uj if Uj > Ui. This makes it possible 
to find the correct ordering of gridpoints (ascending in U) at run-time, ef- 
fectively de-coupling the system of discretized equations. This idea, the 
basis of Dijkstra's classic method for shortest paths on graphs [14] . yields 
the computational complexity of 0(M log M). Tsitsiklis has introduced the 
first Dijkstra-like algorithm for semi-Lagrangian discretization on a Carte- 
sian grid in [40]. A Dijkstra-like Fast Marching Method was introduced by 
Sethian for Eulerian upwind finite-difference discretizations of isotropic front 
propagation problems [29] . The method was later extended by Sethian and 
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collaborators to higher-order accurate discretizations on grids and unstruc- 
tured meshes, in R n and on manifolds, and to related non-Eikonal PDEs; 
see |30j . [33j . and references therein. For anisotropic control problems, the 
discretization (jlOp generally is not causal, but the computational stencil 
can be expanded dynamically to regain the causality. This is the basis of 
Ordered Upwind Methods [34], [35J HI] , whose computational complexity is 
O(TMlogM), where T measures the amount of anisotropy present in the 
problem. 

The idea of alternating the directions of Gauss-Seidel sweeps to "guess" 
the correct mesh point ordering was employed by Boue and Dupuis to speed 
up the convergence in [7]. For Eulerian discretizations of HJB PDEs, the 
same technique was later used as a basis of Fast Sweeping Methods by Tsai, 
Cheng, Osher, and Zhao [39, 43J. Even though a finite number of Gauss- 
Seidel sweeps is required in the Eikonal case, resulting in a computational 
cost of O(M), the actual number of sweeps is proportional to the maximum 
number of times an optimal trajectory may be changing direction from quad- 
rant to quadrant. A computational comparison of fast marching and fast 
sweeping approaches to Eikonal PDEs can be found in [20^ [19] . 

Section 3. Multi-criterion optimal control. 

Unlike the classical case considered in section [2j in realistic applications 
there is often more than one criterion for optimality. For a system controlled 
in C R n , there may be a number of different running costs Kq, . . . ,K r 
and a number of different terminal costs qo, . . . ,q r (all of them satisfying 
assumptions (AO- A3)) resulting in (r + 1) different definitions Jo,---,Jr 
of the total cost for a particular control. A common practical problem is 
to find a control a*(-) minimizing Jq{x, a*(-)) but subject to constraints 
Ji(x, a*(-)) < Bi for all % = 1, . . . , r. 

We will refer to a control minimizing j7o without any regard to constraints 
as primary- optimal. A control minimizing some Jj (for j > 0) without any 
regard to other constraints will be called j -optimal (or simply secondary- 
optimal, when j is clear from the context). A control minimizing J§ subject 
to all of the above constraints on J^'s will be called constrained- optimal. 

For a fixed x E f2, we will say that a control ai(-) is j-dominated by a con- 
trol a 2 (-) if Ji{x,a 2 {-)) < Ji.(x, a%(-)) for all i = 0, . . . , r and Jj{x, a 2 (-)) < 
Jj{x, ai(-)). We will also say that ai(-) is dominated by a 2 (-) if it is j- 
dominated for some j G {0, ...,r}. E.g., the constrained-optimal control 
a*(-) described above might be dominated, but will not be 0-dominated by 
any control. Any control a**(-) dominating a*(-) will have the same pri- 
mary cost Jq and will also satisfy the same constraints; moreover, it will 
even spend less in at least one of the secondary costs J\ , . . . , J r . 

We will say that a(-) is Pareto- optimal for x, if it is not dominated by 
any other control. In that case, its vector of costs corresponds to a point on 
a Pareto-front V(x) in Si JTq; ■ • • ; *Jr space; sec Figure [2lA. 
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Our goal is to find an efficient numerical method for approximating the 
costs associated with Pareto-optimal controls. The efficiency requires solving 
this problem for all x G Q simultaneously. 

We begin by showing how to compute the total cost Ji incurred by using 
a control optimal with respect to a different running cost (section I3.1j> . We 
then describe a recent method due to Mitchell and Sastry for recovering a 
convex portion of the Pareto front by scalarization (section I3.2p . Finally, 
in section 13.31 we describe a new method for solving fully this problem by 
augmenting the system state to include the "budget remaining" in each 
secondary cost. 

Subsection 3.1. Total cost along "otherwise optimal" trajectories. 

We will use Ui(x) to denote the value function with respect to Ji if all other 
costs are ignored. As explained in section [2TTT Ui(x) can be recovered as a 
viscosity solution of the PDE ([5]), if we set K = Ki and q = q%- 

Given a different value function u derived for some other cost J ', we will 
define a restricted set of ^-optimal controls 

A u , x = {o(-) G A | J(x, o(-)) = u(x)} . 

As explained in section \2.1\ if J satisfies the assumptions (AO- A3) , then 
A u x will be non-empty for every x G Q and will contain a single element 
for every x G Q\T. 

We will use Vi(x) to denote a ^7-optimality-restricted value function with 
respect to Jf. 

(11) Vi(x) = inf Ji(aj,a(-)). 

This notation relies on a fixed choice of J and u, and we will specify J 
explicitly in each case to avoid ambiguity. For the cases when J = Jj and 
u = Uj, we will use Vij instead of V{. According to this definition, we also 
have vu = Uj. 

The following properties of Pareto fronts follow from the above definitions. 
The proofs are simple and we omit them for the sake of brevity. 

Property 3.1. Vi(x) > Ui(x) for Vx G Cl and for any choice of J satisfying 
assumptions (AO- A3). 

Property 3.2. Suppose u(x) is the value function corresponding to some 
J that satisfies assumptions (A0-A3). If Ki and qi also satisfy (A0-A2), 
then Vi{x) is lower semi- continuous on Vt and continuous wherever u(x) is 
differentiate. 

Ji > Uj(x), j = 0, . . . ,r} . 



,r. 



Property 3.3. LetU{x) = {(J , . . . , J r ) G R r+ \ 

Then for \/x G £1, 

(1) V{x) C U{x) and 

(2) (vi (x), v ir (x)) G V(x) for all i = 0, . . . 
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FIGURE 1. Cost along "otherwise optimal" trajectories. A simple 
one- dimensional example: Cl = [0, 1], / = Kq = K\ — 1, and qo(0) = 
qi(0) = = 0, but go(l) = 0.4. Similar discontinuities may occur 

even if qo = qi (when Ko K\). 



Property 3.4. If r = 1, thenV(x) C [foo(a;), foi(*)] X [vn{x), v\o(x)]. 



Suppose u is a smooth solution to the PDE ([5]). Then T = and for 
every x € ft, there exists a unique optimal/minimizing control value a* = 
a*(x, V«(i)) € A. The local rate of increase of Ji along the corresponding 
trajectory is Ki(x, a*). This yields a system of (r + 1) linear PDEs 

Vvi(x) ■ f(x, a*) = —Ki(x,a*), \/x £ S7; 
(12) fj(ic) = Va; £ dft; i = 0,...,r. 

This system is coupled to a nonlinear equation ([5]), since a* is generally 
not available a priori unless Vw is already known. 

In the Eikonal case (when A = 5 n _i, f(x,a) = f{x)a and K(x,a) = 
K(x)), the optimal direction of motion is 

Vuix) „ . . /(a?) 

a * = - ; r = -Vu(x) 



\Vu(x)\ y ' K{x) 
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and the system (I12p can be rewritten as 

Vvi(x)-Vu(x) = Ki(x,a*)K(x)/f(x), Vx E Q; 

(13) Vi(x) = qi(x) Vx € d£l; i = 0, ...,r. 

If u is not smooth, the functions Vi may become discontinuous (see Figure 
[1]) and a generalized solution is needed to define i>i(x) at any points x £T. 
Luckily, Bellman's optimality principle provides an alternative definition to 
resolve this ambiguity: 

Vi(x) = lim min {tK(x, a*) + V{ (x + rf(x, a*)) } , 

where A U)X C A is the set of minimizing control values in (j5J). If x £ V, 
then the A u<x consists of a single element and this formulation is equivalent 
to (]12p . Whether or not a* is unique, the above formula yields the lower 
semi-continuity of It can also be used (with a fixed small r > 0) to derive 
a first-order semi-Lagrangian discretization of (|12p . 

We note that the key technical idea employed above (solving a linear 
equation along the characteristics of another PDE) is well-known and useful 
in many applications. A common version of it arises whenever there is a 
need to "propagate a constant" from the boundary along the characteristics 
of some PDE. For example, this is the essential idea behind the "velocity 
extension equation" in [31] and the "escape equations" in [32]. A slightly 
less general version of our equation (|13p (for isotropic costs/dynamics with 
/ = 1 and qi = 0) has also been previously used in |25j . 

Subsection 3.2. Weighted sums method and Pareto fronts. Scalar- 
ization is a popular approach where a point on Pareto-front is recovered 
by minimizing an "aggregate objective function" . The method of weighted 
sums defines that aggregate function as a convex linear combination of the 
original objectives |24j . 

A recent method based on this approach was introduced by Mitchell & 
Sastry for multiobjective optimal control in the case when / = 1 and all 
costs are isotropic |25| . Here we describe a slightly generalized version of 
their method. 

r 

Let A = {A = (Ac • • • , A r ) | = 1 anc ^ ^« — f° r an an d suppose 

A is some mesh discretizing A. 
Given AG A, define 

r r 

K\(x, a) = ^i K i( x , a ); Qx(x, a) = ^2 a )- 

i=0 i=0 

Solve equation ([5]) for K = K\ and q = q\. 

Having found u, solve the system (fT2|) to obtain Vi's. 

The resulting point (vq(x), . . . , v r (x)) will belong to the Pareto front V(x). 

The above procedure is applied repeatedly for all meshpoints A G A to 
obtain a mesh approximating V(x). Since finding each point on the Pareto 
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FIGURE 2. Pareto Front (left) and its reconstruction using the 
"weighted sums" method (right). An optimum found for each A 6 A 
yields a point on a convex part of Pareto front and the vector A will be 
orthogonal to a support hyperplane to the front at that point (shown 
by a dashed line). In weighted sums method, an envelope of all sup- 
port hyperlplanes is used to approximate the Pareto front, but misses 
all non-convex parts of the front |13j . 



front involves solving one non-linear and (r + 1) linear PDEs, any restriction 
of the computational domain f2 C R n is worthwhile. This can be often 
accomplished by finding u±, . . . ,u r first and excluding all x for which Ui(x) > 
Bi for some i G {1, . . . , r}. 

Aside from the computational cost of the above procedure, this approach 
suffers from two usual problems associated with the method of weighted 
sums. First, a uniform mesh on A often results in a highly non-uniform 
mesh on V(x). Secondly, the weighted sum method can approximate only 
a convex part of the Pareto front [13]; see Figure [2)3. This may result in 
selecting suboptimal trajectories. Mitchell & Sastry acknowledge that, for 
non-convex fronts, "this method may fail to detect a feasible path even if 
one exists" , but they report that "non-convexity has not been a problem" in 
their numerical experiments. They further suggest that "neighboring values 
of A can be used to bound the error in the convex approximation" of non- 
convex portions of the front. We believe that the latter procedure can be 
very inaccurate, especially since the Pareto front is frequently discontinuous 
for multiobjective optimal control problems. See Figure [TOl for an example 
of non-convexities actually present in the test-problem introduced in [25J. 

In addition, we note that recovering the entire Pareto front for each x £ f2 
is excessive and unnecessary when the real goal is to solve the problem for a 
fixed list of constraints B\, . . . , B r . If instead of using a mesh A, we attempt 
changing A adaptively, it is not generally clear in what direction should A 
be perturbed in order to satisfy the constraints. 
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Despite the above limitations, this technique can be useful in many ap- 
plications, whenever there is a need to produce at least some of the Pareto- 
optimal trajectories efficiently. E.g., a similar scalarization approach has 
been independently used by Kim and Hespanha in path planning (to mini- 
mize the risk of detection/interception) for groups of UAVs 



Subsection 3.3. An augmented PDE on an expanded state space. 

We propose an alternative approach based on augmenting the system-state 
space to reflect the budget remaining in each of the secondary costs. This 
is a generalization of the idea classically used to recast a Bolza problem as 
a problem with zero running cost by adding an extra dimension to the state 
space [8]. 

Suppose hi G [0, .Bj] is the "budget" remaining in the secondary cost Ji. 
We define an extended state variable x = (x, b%, . . . , b r ) and the extended 
state space Vt e = O X (0, B{) X . . . X (0, B r ). The outflow boundary and the 
inflow boundary of this domain are 

-Bout = {x £ dCl e | x G O, and bi G (0, Bi] for i = 1,. . . , r} and B- in = <9f2 e \£> ou t- 

We note that B OVL t is the part of the boundary where at least one of the 
budgets is at the maximal level (3j such that bj = Bj). For the case r = 1, 
B m coincides with the so called parabolic boundary [15] ; see Figure El We 
will also refer to a feasible subset of B m : 

Bf = {x G B m I x G <90, and bi > qi(x) for i = 1, . . . , r} . 

The extended state at the time t will be denoted as y{t) = (y(t), fti(t), . . . , f3 r {t)) G 
fj e , and the extended dynamics is now defined by 

y'(t) = f(y(t),a(t)); 

(14) fi(t) = -Ki(y(t),a(t)), i = l,...,r; 

(15) y(0) = x = (x,b 1} ...,b r ) Gfi e - 

As before, the exit time corresponding to a particular control is defined 
by ©> but we will use T = a to simplify the notation. 
The total cost of this control is defined as 

(16) J(x,a(-)) = llo K o(y(t)^(t))dt + q (y(T)) if y(T) G B h 

I +oo otherwise. 

The (possibly infinite) value function of the new control problem is defined, 
as usual: w(x) = w(x, bi, . . . , b r ) = inf j(x, a(-)). 

If the value function is smooth, the standard argument (based on Taylor- 
expanding w (y(r)) in the optimality principle), shows that w satisfies the 
PDE 

(17) mmS^K (x,a) + V x w ■ f(x,a) - J2 K i( x ' a )|^| = °' 
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with the boundary conditions 



(18) 



w(x) 



q (x) on Br, 
+00 on B- m \B{ 



However, w can be not only non-smooth, but also discontinuous inside O e 
since the boundary data is discontinuous and no local controllability (in the 
directions b\,...,b r ) is present; see Remark 12. 11 Nevertheless, w can still be 
interpreted as a unique discontinuous viscosity solution [37] of equation (fT71) , 
despite the fact that Soner's "inward pointing" condition is clearly violated 
on B out [26]. 

As in the single-objective case, if the minimization in a can be performed 
analytically, the augmented PDE can be rewritten in a simpler form. For 
example, in the fully isotropic case, 



(19) 



dw 



K (x) - \V x w\f(x) - Y,Ki(x)— = 



i=l 



which in the case r = reduces to the usual Eikonal equation ([7]). 



Bi 



bi < ui(x) 



w(x, 61) = uo(x) 



w(x, 61) = v ,i{x) 







FIGURE 3. A sketch of the domain f2 e for the case r = 1. The thick 
line shows the inflow boundary Bin ■ The dashed line is the graph of ui (x) 
(the minimum 61 needed to leave fi starting from x). Bt is the portion 
of Bin above the dashed line. The dotted line shows the points where 
bi = vifi(x). At that level and higher, the primary-optimal controls are 
feasible and w(x,bi) — uo(x). Above the dotted line, the constraint is 
slack. The PDE (|17[) has to be solved numerically in the shaded region. 



16 



AJEET KUMAR AND ALEXANDER VLADIMIRSKY 



The following properties of w follow from the above definitions. The 
proofs are simple or standard and we omit most of them (except for the last 
three) for the sake of brevity. 

Property 3.5. The value function w : Cl e i— > (i£U {+00}) is lower- semicontinuous; 
see [37] and references therein. 

Property 3.6. The value function w is monotone non-increasing in each 
hi- As a result, ifb > c componentwise, then w(x,b) < w(x,c). 

Property 3.7. If for some i, bi < mix), then w(x,b) = +00. 

Property 3.8. w(x, b) > uq(x) = w (x, v\q(x), . . . , v r o(x)) . 

Property 3.9. w is Lipschitz continuous along every characteristic. 

Property 3.10. The characteristics of PDE (JT7J have the following prop- 
erties: 

(1) All characteristics are monotone-increasing in all bi's. 

(2) Projections of characteristics on £1 yield constrained- optimal trajec- 
tories. 

(3) Any characteristic starting on Bt leads into Q e . 

(4) Any characteristic starting on B ou t leads out of Q e . 

We will say that a control a(-) is feasible for x = (x, b) if 3%{x, a(-)) < bi 
for all i = 1, . . . , r. We will say that x € Cl e is a feasible point if it has at 
least one feasible control (i.e., if w(x) < 00). We will say that a point x is 
i -tightly- constrained if for every constrained-optimal control «*(•) we have 
Ji(x,a*{-)) = b{. Otherwise, we will call x i-slack. We will also say that x 
is totally slack if there exists a constrained-optimal control a*(-) such that 
Ji(x, o*(-)) < bi for alH = 1, . . . , r. 

Property 3.11. The point (w(x,b),b) lies on the Pareto front V{x) in 
R r+1 if and only if (x, b) is i-tightly- constrained for all i = 1, . . . , r. 
Moreover, if (bo,...,b r ) € V(x) and bi < Bi for all i = l,...,r, then 
w(x,bi,. . . ,b r ) = b . 

Proof. Suppose a*(-) is the constrained-optimal control for x = (x,b). If 
x is i-slack, then Ji{x,a*{-)) < bi and this does not contribute a point on 
Pareto front. 

If x is i-tightly-constrained for all i = l,...,r, then <!*(•) cannot be 
dominated by any control. (If it were O-dominated, that would contradict 
its constrained-optimality. If it were i-dominated for % > 0, that would 
contradict the fact that x is i-tightly-constrained.) 

Finally, suppose a(-) is a control that realizes the cost vector (bo, . . . , b r ) £ 
V(x) with bi < Bi for all « = 1, . . . ,r. By the definition of w, we have 
w(x, 61, . . . , b r ) < Jo{x,a(-)) = 6o- If we had w(x, b\, . . . , b r ) < bo this 
would imply the existence of a control 0-dominating a(-), which would con- 
tradict its Pareto-optimality. So, w[x, b\, . . . , b r ) = bo. □ 
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Property 3.12. If the point (sci, 6) is totally slack, then w is locally Lipschitz- 
continuous in its first argument. 

Proof. First, we note that there exists some open neighborhood of x\ such 
that (x2,b) is totally slack for every X2 from that neighborhood. Otherwise 
there would exist a sequence of i-tightly-constrained Xj's converging to x%, 
which can be used to show that (x\, b) is i-tightly-constrained as well. 

If (sci, b) is totally slack, then, starting from any point X2 close to Xi, the 
i-th cost of reaching x\ is bounded above by |cci — £C2 1 m ax(iQ/|/|), where the 
maximum is taken over all x and all a such that f(x, a)/Ki(x, a) £ dVi(x). 
By the assumption (A3), V{ contains the origin in its interior; so, there exists 
a constant L such that the i— th cost of reaching x\ from X2 is bounded by 
L\x\ — x<i\ for all i. Suppose we travel from (x2, b) to (xx, b), where b <b. 
If X2 is sufficiently close to xi, then (xi,fo) is still totally slack, and any 
constrained-optimal control for (cci, b) will be still feasible for (xi, 6). Thus, 
w(x2,b) < L\x\ — xi\ + w(xi,b). To complete the proof, we repeat the 
argument switching the roles of X\ and X2- □ 

Property 3.13. If x = (x,b) is totally slack, £»*(■) is a constrained- optimal 
control for x and y*(t) is the corresponding trajectory in f2, then y*(t) is 
also a characteristic of problem © with K = Kq and q = qo- 

Proof. Briefly, if (x, b) is totally slack, then any sufficiently small pertur- 
bation of a*(-) will also be feasible. Since a*(-) is constrained-optimal, 
the function y*(t) is a local minimizer of j7o and will be a solution of the 
Euler-Lagrange equation in R n (see, e.g., p2]). By Pontryagin's maximum 
principle, it will also be a characteristic of the corresponding HJ PDE on 
Cl. This is an interesting fact, since the characteristics of ([5]) yield locally 
primary-optimal (unconstrained) trajectories. □ 

Property l3.1Q[ l is the basis for explicitly causal discretizations of the aug- 
mented PDE, which enables an efficient numerical treatment (by "marching" 
in any direction hi). Properties 13.71 and 13.81 can be used to reduce the com- 
putational domain, as shown in Figure [3] and further explained in section 
13.51 Property 13.111 can be used to extract the relevant part of the Pareto 
front from the values of w. 

Subsection 3.4. Numerical method for the augmented PDE. We 

consider a Cartesian grid X on Cl e . For simplicity, we will assume the same 
grid spacing h in all spatial dimensions and grid spacing A&i, . . . , Ab r in each 
of the constraint /secondary cost directions. Let h = max{/i, maxj A6j}. Our 
goal is to construct an approximate solution W converging to the lower semi- 
continuous value function w as h — > 0. 

If the minimization in a can be performed analytically (e.g., in the fully 
isotropic case of equation (|19h ). a natural Eulerian framework scheme may 
be obtained by using upwind finite differences to approximate the partial 
derivatives of w. However, this approach, even when feasible, would be 
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subject to CFL-type stability conditions, which would potentially have a 
significant impact on the computational cost of the scheme. (A simple ex- 
ample to illustrate this: suppose the / and Kq are isotropic, r = 1, q\ = 
on d£l, and K\ = 1, making J\ be the total time to <9f2 along a given trajec- 
tory. If upwind finite differences are used in a given &i-slice to approximate 
V x w and the forward-difference is used to approximate this results in 
a scheme suitable for explicit forward "time" -marching in &i-direction, but 
not surprisingly requires a standard CFL condition maxx f(x) < h/Ab\ for 
stability.) 

Instead, we chose to implement a semi-Lagrangian discretization of the 
augmented PDE (|17p . In addition to improved stability properties, the 
resulting scheme is also easy to extend to unstructured meshes in £l e and 
is applicable when the minimization in a cannot be handled analytically 
(which is often the case for control-theoretic problems). Our discretization 
is a variant of ([9]). Given a point x = (x, b) = (x\, . . . , x n , bi, . . . , b r ), we 
define a new system state (x, b) as follows: 

(20) x = x + T a f(x,a); % = bi - r a Ki(x , a) , for i = 1, . . . ,r. 

Here the obvious dependence of the new state on the choice of control a is 
omitted for the sake of notational simplicity. 

(21) 

W(x, bi, . . . , b T ) = min ^T a Ko(x, a) + VF(a:,6)j, V gridpoints (x, b\, . . . , b r ) ^ B; 

To ensure that the discretized equations allow efficient marching in the 
direction b\, it is sufficient to take 

(22) r a = Ah/K^x,^, 

which guarantees that b\ = b\ — Ab\ for any choice of a. This means that 
the new position (x, b) will be in the previous "&i-slice", in which the values 
of W were already computed. 

Of course, (x,b) will usually not be a gridpoint and W(x,b) has to be 
approximated using the values of W at nearby gridpoints. Our implemen- 
tation uses a standard tensor product of linear interpolations in all x and b 
variables; see, e.g., [9]. 

For example, if n = 2 and r = 1, this yields a bilinear interpolation. Sup- 
pose x lies in a cell with vertices x\, . . . , x^ (enumerated clockwise, starting 
from the lower left corner). Let (71,72) = (x — X\) /h. Then 

W(xM) = 71 (72^x3, 61) + (l-72)W(x 2 ,h)) 
23 V , „ ' ~ \ 

+ (1 - 71 ) (72^(^4, 61) + (1 - l2)W{x x M)) ■ 

Remark 3.14. If the value function is smooth on the cell, the resulting 
interpolation error is 0(h 2 ), where h is the size of that cell (i.e., h = 
max{/i, maxAfrj}). However, in the worst case, w may be discontinuous, 
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resulting in a 0(1) interpolation error. (We note that the property 13.91 
ensures the Lipschitz continuity of w on the characteristic itself, but the 
interpolation cell containing x might still include a discontinuity.) 

The convergence of semi-Lagrangian schemes to discontinuous viscosity 
solutions is more subtle. In [U [5] Bardi, Falcone and Soravia have proven 
that, on any compact subset where w is continuous, the semi-Lagrangian 
approximation converges to w uniformly, provided h/r — > as r — > 0. 
The resulting schemes were successfully used to approximate discontinuous 
viscosity solutions both in the context of infinite-horizon optimal control 
and in differential games. For viscosity solutions continuous on the entire 
domain, a different proof [5] yields the convergence (and the convergence 
rate estimates) under a milder assumption that h/r remains bounded as 
r -> 0. 

Since we are setting r for each x and a separately, the above condition 
for convergence to discontinuous solutions becomes: fo/(minr a ) — > as 
(maxT a ) — > 0. If r a is selected according to (|22j) . this condition requires 
that the grid refinement should be performed in such a way that h = o(A£>i). 
(An alternative approach would be to keep h = O(Abi), but pick r a so that 
T a K±(x, a) = b\ — b\ = mA&i, where the integer m — > oo and mAbi — > 
as h — > 0.) However, for the examples considered in this paper, the 
numerical evidence suggests convergence even with m = 1 (i.e., as prescribed 
by formula (|22p ) and with h = O(Abi). The discontinuities are smeared in a 
narrow band with the width of that smearing proportional to Ab±; e.g., see 
the convergence study for a simple example in section 14.11 We note that any 
higher order accurate semi-Lagrangian scheme would need to use a larger 
stencil to interpolate W(x, b), and an ENO or WENO reconstruction would 
be needed to handle the discontinuities. 

Even though it is not necessary in principle, our current implementation 
assumes the geometric dynamics (defined in section 12. ip . The minimization 
in formula (|2ip is performed numerically using the standard "golden section 
search" algorithm. Once W is computed, optimal controls and trajectories 
are recovered by following the characteristics of PDE (|17p numerically. 

Subsection 3.5. Reducing the computational domain & initializa- 
tion. Given the high dimensionality of Cl e , any reduction of the domain 
is likely to result in substantial savings in the computational cost. Two 
such reductions are obviously worthwhile; see Figure El In both cases we 
recover a surface consisting of special characteristics of (I17p by efficiently 
solving other PDEs on lower-dimensional domains. We note that a similar 
approach was previously proposed in [42] for time-optimal control in the 
case of time-dependent dynamics. 

First, by property 13.81 if w(x,b) = uq(x), a primary-optimal control is 
already feasible and further increase in secondary budgets will not yield any 
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reduction of w. The n-dimensional surface on which this happens can be 
found a priori (by numerically approximating functions vm for i = 1, . . . r). 

Secondly, only a subset of fl e is feasible: if the initial budget-vector b is 
insufficient to reach d£l starting from x, then w(x, b) = +00. We formally 
define the Minimal Feasible Level (MFL) as a graph of the function 

w\(x, 6 2 , ... , b r ) = min{6i | w(x, b\, b 2 , . ■ . , b r ) < +00} . 

As described below, the MFL can be efficiently recovered for any r by solving 
a sequence of PDEs on lower-dimensional domains. However, for r = 1, this 
task is particularly simple, since in that case the MFL coincides with the 
graph of u\(x), and the latter can be approximated by solving the PDE ([5]) 
numerically on Q. As explained in |37j . the augmented PDE (|17p then has 
to be solved on the epigraph of u\. Thus, using the results in section [3TTT 
the valu^l of w on the MFL is provided by Vq\(x). 

Remark 3.15. To represent MFL on the grid X, our current implementa- 
tion uses the smallest integer i such that u\(x) < iAbi and then initializes 
W(x,iAb\) = vqi(x). This procedure is conservative (in the sense that 
it overestimates the minimal necessary budget), but it also introduces ad- 
ditional errors of order 0(Ab\); see the convergence study in section 14.11 
A better implementation can be built by locally altering the interpolation 
procedure near the MFL. 

For r > 1 the approximation of the MFL is more subtle. In [37J, So- 
ravia suggested solving an augmented PDE on the set {(x, b) £ Cl e | 6j > 
Ui(x), i = l,... r}. However, we note that w need not be finite everywhere 
on that set. This is already evident for r = 2. Indeed, a control optimal 
according to J\ usually incurs a higher J^-cost than the control optimal 
according to the latter (and vice versa). As a result, w (x, u%(x), U2(x)) 
is usually +00. This situation is illustrated in Figure H] (i.e., the region 
between the dashed and solid lines is not feasible). If 62 > V2i(x), then 
w(x,ui(x),b2) = Vqi(x); this allows to recover a part of the MFL quite 
easily (the thick solid lines in Figure Hj). However, recovering the rest of the 
MFL requires more computational effort. 

Suppose r = 2 and W\(x, 62) is the minimum budget b\ sufficient to reach 
dVi from x with J2 < 62. It is clear that the PDE for w should be solved 
on the epigraph of w\. Since w\ optimizes J\ subject to a constraint on 
J2, we note that w\ plays the same role for the pair of costs (J7i, ^2) as w 
played for the triple (Jo, J\, ^2)- Thus, the same argument used in section 
13.31 shows that w\ can be recovered as a viscosity solution of a PDE similar 



We emphasize that the MFL does not really provide any additional boundary condi- 
tions for the equation (I17|l (since this surface consists of characteristics of that PDE). But 
for a semi-Lagrangian discretization, the numerical values on the MFL are needed - for 
an x just above that surface, the optimal x might well lie in a cell one of whose corners is 
on the MFL. We note that any x falling below the MFL is obviously non-optimal; when 
such situation arises during the minimization in equation (|21[) . we use W(x) — +00. 
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FIGURE 4. Two secondary costs diagram for a single point x £ Q. 
The boundary of the feasible set (x's MFL) is shown by a solid line. To 
solve the PDE using an efficient semi-Lagrangian scheme, it is necessary 
to find the entire MFL (including its curved and possibly non-convex 
part) and pre-initialize w on that boundary. The dotted line shows the 
boundary of the set where the unconstrained (primary-optimal) controls 
are feasible. 



to (fT7|) . but posed on SI x [0,-62]; i-e., 



(24) 



mjn <j Ki(x, a) + V x Wi • f(x, a) 



K 2 (x,a 



0. 



To solve the latter PDE efficiently on an (n + l)-dimensional domain, we 
first need to approximate its own n-dimensional "MFLi". Luckily, this 
is equivalent to the one-secondary-cost problem already considered above. 
Thus, this new MFL\ is the graph of u 2 {x) and the value of w\ on the 
MFL\ is provided by v\2{x). While computing w\, we can also integrate 
Kq along wi's optimal trajectories (as described in section [37T]) to find the 
values of w on its MFL (i.e., on the graph of w\ in Ct x [0, Si] x [0, B2]). 
This is the approach used to compute the example in section 14.81 

For the general case (r > 2), the same procedure can be applied re- 
cursively. Let uii(x, foi+i, . . . , b r ) be the minimum budget bi sufficient to 
reach d£l from x with Jj < bj for all j > i. According to this definition, 
w r (x) = u r (x) and, for \fi < r, W{ is finite only on the epigraph of Wi + \. 
Focusing only on the costs Ji, . . . J T , the argument of section [3731 shows that 
Wi is the (unique lower semi-continuous) viscosity solution of 



(25) 



mm ^ Ki(x, a) + V x Wi ■ f(x, a) 



j=i+i 



dwj 
db 3 
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on O x [0, .Bj-i-i] x ... x [0,S r ]. In this (n + r — z)-dimensional domain, 
the MFL for W{ (i.e., the MFLi) is provided by a graph of Wi+i, and the 
values of Wi on MFLi can be computed by integrating K{ along the optimal 
trajectories of itfi+i. Finally, the graph of w\ is the MFL for w = Wq, and 
the latter solves the PDE (1T7|) . 

Remark 3.16. Both techniques presented in this section allow a significant 
reduction of the computational domain. The net effect of the first technique 
(based on considering only those gridpoints where w(x, b) > uq(x)) depends 
on how large the budget bounds B{S are compared to the values of Vio(x). 
If r = 1, the percent of excluded gridpoints can be found by averaging the 
ratio [max{£>i — v\o(x), 0} / B{\ over ft. 

The net effect of the second technique (based on finding the MFL) depends 
on the degree of "primary-non-optimality" of "secondary optimal" trajecto- 
ries. If r = 1, the percent of gridpoints additionally excluded (after the MFL 
is computed) can be found by averaging the ratio [min {Bi,u\(x)} / min {Bi,v\o( 
over f2. Based on the experimental evidence, the resulting efficiency gains are 
quite substantial. E.g., in the examples of sections l4.1l and l4.2.1l the MFL 
allowed additional exclusion of 50% and 93% of the remaining gridpoints 
respectively. In the example of section [4~8| where r = 2, the corresponding 
reduction was 76%, confirming that the recursive procedure to compute the 
MFL when r > 1 is worthwhile. 

Subsection 3.6. Selecting r & optimal ordering. Formula (|20p for the 

new state after using control a for r a seconds is based on the assumption 
that, for every fixed a, / and all K^s are approximately constant near x. 
There are two advantages to this formula. First, it allows to compute (x, b) 
very quickly. Second, taking r a = Ab\/Ki(x, a) ensures that the new state 
is in the previous &i-slice; i.e., b\ = b\ — Ab\. This allows for the explicit 
causality and marching in the b\ direction. In section 14.91 we wn l refer to 
this implementation as Algorithm 1 (see Figure [5A) . 

There are also two obvious drawbacks to this discretization. First, this 
linear approximation is poor wherever / and i*Q's vary significantly near x. 
Second, the local error in formula (|2ip is generally 0(t%) even for a smooth 
w and perfectly known W(x,b); thus, it would be preferable to select the 
smallest r a that still allows explicit marching in the direction b\. 

To address the first of these drawbacks, we have implemented Algorithm 2. 
Given the current state x = (x, b) and a particular control value a, we start 
with the ray from x in the direction [f(x, a), —Ki(x, a), . . . , —K r (x, a)]. We 
find the first intersection of that ray with an (n+r— l)-dimensional cell of the 
grid. If that cell is fully in the previous &i-slice, we can interpolate, as in Al- 
gorithm 1. Otherwise, we re-evaluate /, and K^s at the new point and follow 
the new ray, repeating the process until we reach the previous &i-slice (see 
Figure 03). This procedure is computationally more expensive than Algo- 
rithm 1, since we might need to traverse through several (n + r)-dimensional 
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h - Ah 
(A) 



x 



X\ 



X-2 



X 3 



X\ 



61 - A61 = b[ 



X 



X\ 



FIGURE 5. Choices of r illustrated for r = 1 and n = 1. Let x = 
(x4, 61). A trajectory corresponding to a particular a is shown by a dot- 
ted line. Algorithm 1 is illustrated on the left: set ra = Abi/Ki(x,a) 
and assume that / and KiS don't change. This ensures bi = b\ and gives 
a direct formula for x 6 [asi, X2] ■ Algorithm 2 is illustrated on the right: 
/ and Ki's are re-evaluated at each intersection with 1-dimensional cells. 
Algorithm 3 is based on the fact that, if W(xz, b\) has been already com- 
puted, then there is no need to continue beyond the first intersection. In 
that much smaller r a already guarantees the causality, x = X3, 

and 61 £ [b[, 61]. 



X 2 



x 3 x 4 



cells before finding x, especially when Ab\/min(K\) > h/min\f\. How- 
ever, it is much more suitable for problems with rapidly varying / and K^s. 
(Indeed, in section H~9l we consider several numerical examples, where these 
functions are actually discontinuous in x.) 

To alleviate the second problem, we have implemented Algorithm 3. The 
idea is to select the smallest r a needed for the causality. Our implementation 
uses the smallest r a sufficient to guarantee that (x, b) lies in an (n + r — 1)- 
dimensional cell, all vertices of which have already had values of W previ- 
ously computed (making interpolation possible). This will certainly be the 
case if that cell lies fully in the previous 61-slice (as in Algorithm 2), but it 
may also happen earlier; see Figure [SB for an example. We note that all 
relevant intersection points are already computed in Algorithm 2. If one of 
them is suitable in the above sense, we don't need to continue beyond it 
(thus reducing both the computational cost and the local truncation error) . 
In the worst case, we still need to trace this piecewise-linear trajectory up 
to its intersection with the previous 61-slice (as in Algorithm 2). As a re- 
sult, the computational cost of this Algorithm is always somewhat less than 
that of Algorithm 2 though obviously higher than that of Algorithm 1; see 
section 14.91 

Remark 3.17. We note that both the accuracy and the efficiency of Algo- 
rithm 3 clearly depend on the order of computing Ws within the same 61- 
slice. (The location of x in Figure will be different depending on whether 
we have already computed W(x^,bi) prior to W(x^,bi).) Our current im- 
plementation uses a simple lexicographic ordering in secondary costs and 
then in spatial directions: assuming that x = (x, b) = (x%, . . . , x n , b%, . . . ,b r ), 
we compute within each 61 -slice in the direction b 2 , within each (61, &2)-slice 
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in the direction 63, ... , within each (6)-slice in the direction xi, within each 
(6, xi)-slice in the direction x%, etc. While the preference for the secondary- 
cost (over spatial) directions is clearly motivated by the explicit causality 
(i.e., Ki S are positive), our fixed ordering of the spatial directions is arbi- 
trary and hardly optimal. In the future, we would like to investigate the 
effect of using different orderings in spatial directions (e.g., in the r = 1 
case, sorting ic's by W values found in the previous 61-slice). 

We emphasize that the goal of finding a good ordering within a 61-slice is 
simply to reduce the local truncation error and to speed up each 61-step. 



We illustrate our numerical method with several examples. Our approach 
requires to choose a primary objective function and treat other objectives 
as secondary. All feasible controls will satisfy the constraints Ji < Bi for 
i = 1, . . . , r. We then minimize Jq along these feasible paths. 

All of the examples in this section are computed for a two-dimensional 
system state and we assume that Cl = [0, 1] x [0, 1]. In all of our examples 
(except the convergence study in section 14. ip , the goal is to reach a single 
target point. We will thus assume that qi = at the target and % = +00 on 
the rest of d£l for i = 0, . . . , r. Wherever not explicitly specified, the target 
is at the point (1, 1). 

Most of the examples (except for section l4.2.2p are fully isotropic in cost & 
dynamics, and the results are obtained by solving a variant of equation (|19() . 
For the isotropic case, the dynamics of the system simplifies as y'(t) = 



and a 6 [0, 2ir] is the control parameter. 

In each test problem, we plot constrained-optimal paths for different re- 
source vectors b. We also show secondary -optimal trajectories (recovered 
from UiS, for i = l,...,r) even in cases when these trajectories do not 
satisfy other integral constraints. 

We note that the observability examples in subsections l4.5l - l4T8l use piece- 
wise continuous running costs and/or speed functions. Even though this 
violates assumption (Al), the value function is also well-defined in this case 
and the semi-Lagrangian discretization (based on discretizing Bellman's op- 
timality principle) appear to converge to it correctly. A detailed discussion 
of viscosity solutions to HJB PDEs with discontinuous Lagrangian can be 
found in [38J. A different class of observability- influenced path planning 
problems is considered in |10| . 

All Figures presented in this section were computed using Algorithm 1 
(in subsections 14. II - 14. 3\ where the running costs and speeds are continuous) 
and Algorithm 3 (in subsections 14.51 - 14. 8|) . 



Section 4. Numerical Experiments. 




where y(0) = x G SI C R 2 and / : O 1— > R is the speed 



Subsection 4.1. A simple example: convergence study. We use a two- 
dimensional generalization of example in Figure Q] to test the convergence 
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of our method. We assume that / = 1, and Kq = K\ = 1. Unlike in all 
other examples of this section, here we assume that there are two possible 
exit points on the boundary Ax = (0,0.5) and A 2 = (1,0.5) and define 



qo(x) 



if x = Ad 
if x = A 2 ; 
otherwise. 



if x = A%] 
if x = A 2 ; 
otherwise. 



As a result, J\ is simply the pathlength of the chosen trajectory, all (con- 
strained optimal) trajectories are straight lines leading to Ax or A 2 , and the 
analytical expression for the discontinuous solution w is readily available: 



w(x, bx 



\x - 
\x 

oo, 



A 2 \, 

Ail + 1.5, 



if | sc — A 2 \ < bx; 

if \x — Ai| < bx < \x 

otherwise. 



A, 



We use this example to test the convergence of our method numerically. 
(See Figure [6j) 



Contour plot of Constrained Value function at slice no: 30 (401 *401 *41 ) 
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Lqo error computed on Q'(Abx) 


~\Abx 
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1/10 


1/20 


1/40 


1/40 


0.0039922 


0.0114681 


0.3558851 


1/80 


0.0016614 
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1/160 


0.0007978 
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0.0038505 
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0.0004434 
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0.0002202 


0.0002217 


0.0001994 



FIGURE 6. Convergence to a discontinuous solution (section 14. ip . 
Left: the level sets of w(x, 0.75). Right: the L\ and Loo errors computed 
for various (h, Abi ) combinations. 

We recall that there are typically four different sources of error in semi- 
Lagrangian schemes: (1) due to approximating the boundary conditions 
(in our case - also the approximation of the MFL) on the grid; (2) due 
to interpolation at the foot of the characteristic (e.g., as in formula ([23|) ~1: 
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(3) due to approximating the characteristic and (4) due to approximate 
integration of the solution along that characteristic. 

Since in this example the characteristics are straight lines and Kq is con- 
stant, the last two of these four sources of error are absent. Of the remaining 
two, the first should clearly be decreasing with A&i (since errors of order 
0(A&i) are introduced when representing the MFL on the grid; see Remark 
I3.15p . Assuming that the solution is smooth, each interpolation error is 
0{h 2 ) (due to our use of bilinear interpolation (f23jl ). However, the num- 
ber of interpolations is inversely proportional to A&i and the cumulative 
effect of interpolation errors is larger when the total number of &i-slices in- 
creases. Thus, if h is held constant, decreasing A6i will eventually result in 
an increase in the overall error; this can be seen in the top two rows of the 
corresponding Li-errors table. Of course, in more general problems, where 
characteristics are not straight lines and -fQ's are not constant, the third 
and fourth sources of errors would normally prevent this phenomenon. 

To verify the convergence, we include two tables of error measurements 
for various (h, A&i) combinations. The first table shows errors measured in 
L\ norm on Qf (the subset of £l e on which w is finite though not necessarily 
continuous). For discontinuous solutions, convergence in Lqq norm is gen- 
erally possible only for numerical methods that explicitly track the location 
of discontinuities. Since no such explicit tracking is performed here, we can 
only demonstrate Loo-convergence away from discontinuities. The theoret- 
ical results in [H [31 [5] guarantee that, if h = o(r) as r = Ab\ — > 0, then 
the semi-Lagrangian schemes uniformly converge to the viscosity solution w 
on any compact subset on which w is continuous. For each value of b\, we 
define the discontinuity set of w in the corresponding &i-slice: 



and for each A&i we define a subset of S7 e to study the convergence: 



The second table shows errors measured on Q'(Abi). Since the width 
of the "excluded band" is defined for each column of the second table sep- 
arately, the initial errors in the second and third columns are significantly 
larger, but quickly decrease as h — > 0. In this example, our numerical results 
suggest the convergence even for h = O(r). 

Subsection 4.2. Fastest paths (with restriction on pathlength). Here 
we consider two different examples (one isotropic with obstacles, the other 
anisotropic without obstacles, both inhomogeneous), in which the goal is 
to minimize the time-to-target subject to constraints on the maximal al- 
lowable pathlength. In the absence of obstacles, we use Kq(x) = 1 and 
Ki(x, a) = \f(x, a) | to ensure that j7o and J\ are respectively the time and 
the pathlength along the corresponding trajectory. If obstacles are present, 



Z\ = {(x,b\) G r2 e | w is discontinuous at (x,b\j) 
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we set Kq = +00 inside them, to ensure that all trajectories passing through 
them have infinite cost J§. 

Subsubsection 4.2.1. Isotropic dynamics/cost in the presence of obstacles. 
Here we suppose that the dynamics is isotropic; i.e., f(x, a) = f(x)a, where 
a G S\ = A. For every x = (x, y) G Q, outside of obstacles, we will assume 
that f(x,y) = 1 + 0.8 sin(47rx) sin(67n/) and Ki(x,y) = f(x,y). Consider 



primary: Navigation time, secondary: Path length 




.25 .50 .75 1.00 



Figure 7. Fastest paths, isotropic dynamics (constrained 
by pathlength) and obstacles. Optimal trajectories and the 
level curves of uo. 

a goal of finding the fastest path to the top right corner (1,1) subject to 
restriction on a maximum allowable pathlength and in the presence of three 
rectangular obstacles. The computations are performed on a 301 x 301 x 
301 grid and several trajectories are shown in Figure [7] superimposed on 
top of the level curves of uq. Secondary-optimal trajectory (dotted line) is 
the shortest path. Bold solid line shows a trajectory computed for a large 
maximum pathlength. The constraint is slack in this case and the resulting 
trajectory is in fact primary-optimal; hence its orthogonality to the level 
curves of uq (since for unconstrained isotropic problems the characteristics 
coincide with the gradient lines of uq). The dashed line shows a constrained- 
optimal trajectory for a smaller budget (binding constraint). 
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Subsubsection 4.2.2. Anisotropic inhomogeneous dynamics. The following 
example of anisotropic dynamics was previously used in |34(. [35] , where it 
was motivated by problems in anisotropic seismic imaging. This belongs to 
a class of anisotropic geometric dynamics problems, i.e., f(x, a) = fix, a)a, 
where a € Si = A. Suppose C : [0, 1] h- ► R is a smooth function. We are 
interested in defining / so that, for every x = (x, y) G 12, the "velocity 
profile" V{x) = {f(x,a)a \ a € Si} is an ellipse whose major/minor semi- 
axis have lengths F 2 and Fi respectively and the major semi-axis is aligned 
with the graph of C (i.e., parallel to its tangent at the point x). This is 
attained by setting 



P 

q 



f(x,a) = F 2 \1 + 
The derivation can be found in 35 and 




-1/2 



where 



P 

q 



!+(£(*))' 



WeuseC(a;) = 0.1225 sin (4nx), 



primary: Shotest time, secondary: Path Length Boundary of "Constrained Domain" 

(shown as a thick line) 




.25 .50 .75 1.0 .25 .25 .50 .75 1.0 



A B 

Figure 8. Anisotropic dynamics, fastest paths (constrained 
by pathlength) . Left: level sets of w (x, 1.5) (corresponding to 
level sets of uq in Figure 6 A of [34]) and optimal trajectories 
from (0.1,0.1) to (0.5,0.5). Right: level sets of w(x, 0.6). 



F% = 0.8, Fi = 0.2 and compute the optimal (fastest) trajectories to the cen- 
ter of O. The unconstrained problem leads to an anisotropic static Hamilton- 
Jacobi PDE for no, which can be efficiently solved using Ordered Upwind 
Methods; the level sets of no for the above parameters can be found in Figure 
6 A of [34J). Here we compute pathlength-constrained min-time trajectories 
for the same example on a 201 x 201 x 301 grid. Figure OA. shows optimal 
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trajectories to the center starting from the point (0.1, 0.1) superimposed on 
the level sets of w(x, 1.5) . This is a large "budget": since b\ = 1.5 > v\q{x) 
for \/x £ f2, we have w(x, 1.5) = uq(x) and all unconstrained-optimal tra- 
jectories are feasible. (We note that these optimal trajectories do not follow 
the gradient lines of uq since the speed f(x,a) is anisotropic.) On the 
other hand, for b\ = 0.6, parts of the domain are not reachable (note that 
w = +oo in the corners of the domain in Figure 03)- The thick line sepa- 
rates the part of O where w(x, 0.6) = uq(x); since the starting point (0.1, 0.1) 
is outside that set, its constrained-optimal trajectory is different from the 
unconstrained-optimal. 



Subsection 4.3. Optimal flight-path: minimizing weather threat. 

Here, we consider a test example introduced by Mitchell and Sastry [25j : 
finding an optimal path for an airplane flying from one city to another. In 
one of their examples, weather threat is the primary objective to be mini- 
mized while the fuel minimization is the secondary objective. The running 
fuel cost K\ and the speed of the airplane / are both taken to be unity in 
this example. Hence, the secondary objective, fuel cost is same as the path 
length. 

As mentioned in [25J , weather threat at a location, intuitively, is the prob- 
ability of encountering a storm. Assuming the probability at all locations 
to be independent, the total probability of encountering a storm during the 
flight is one minus the product of the probabilities of not encountering the 
storm at all locations along the path. As we require running costs of integral 
type, the logarithm of these probabilities are used to make weather-threat 
cost additive in nature. The numerical parameters for this simulation are 
listed in table 1. Figure [9] shows optimal paths corresponding to different 
bounds on fuel budget, as well as the contours of weather threat cost. The 
latter is taken to be unity everywhere in Cl except within the two rectangular 
bars where it jumps to a relatively higher value discontinuously. The mag- 
nitude of weather threat is 12 in the darker part of each of the rectangular 
bars and 4 in the other part of the bar. Weather threat function is further 
smoothened to remove the discontinuity and make it Lipschitz continuous. 

As expected, the secondary-optimal path (the shortest path) follows a 
straight line between the two cities. Bounds on fuel budget associated with 
other constrained-optimal paths are shown in the legend. The Figure shows 
that, as the fuel budget increases, the airplane chooses a path through the 
region less susceptible to weather threat. These optimal paths match closely 
those in [25j . 

Figure [10] shows the Pareto front for this example. Pareto front is gen- 
erated using two different approaches. In the first approach, fuel budget 
is taken as secondary objective while in the other approach weather threat 
is the secondary objective. As h and Ab\ decrease, these two versions of 
Pareto front look more and more similar. The non-convex parts of this front 
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have not been found in [25]; this illustrates the advantage of our approach 
as compared to the "weighted-sum scalarization" based techniques. 



primary: weather threat, secondary: fuel budget 
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Figure 9. Path planning for an airplane striving to be away 
from the region having high weather threat probability. The 
optimal paths are shown for the same constraint values used 
in Figure 3 of [25] . 



Grid points in system-domain 


201 x 201 


Grid points along secondary budget 


401 


Speed of the vehicle 


1 


Primary cost (Weather threat) 


explained in 4.1 


Secondary cost (Fuel rate) 


1 


Starting point for optimal path 


(0.1, 0.1) 


Target point for optimal path 


(0.9, 0.9) 



Table 1. Numerical parameters for the optimal flight-path 
example (Figured]). 



Subsection 4.4. Computation of non-visible region. The remaining 
examples deal with robotic navigation in domains with obstacles in the 
presence of friendly and adversarial observers. The observer's position is 
assumed to be known and static and O is split into parts directly visible 
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Pareto Front : 331 by 331 grid size 




2 1 1 1 1 1 1 1 1 l_ 

1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 

tuel budget 

Figure 10. Pareto front for a particular destination for 
the airplane problem. A convexification of this front can be 
found in Figure 4 in |25j . 

(f2„) or invisible (Oj) to each observer due to obstacles. The observabil- 
ity running cost is then set to be high on Q v and low on f2j for an enemy 
observer (and the opposite of this for a friendly observer). 

To enable the computation of constrained-optimal paths, we need to first 
determine £l v and fij. While this can be accomplished by many methods, we 
use the efficient technique based on a Fast Marching Method and described 
in [30]. 

We first solve the Eikonal equation \Vipi(x)\f(x) = 1 with / = 1 on Q 
and the boundary condition ipi = at the observer's location. We then solve 
the Eikonal PDE for ip2 with the same boundary condition but with / = 
inside the obstacles. The region with ip2 > ipi defines the non-visible region. 
In practice, we use a (heuristically adjusted) threshold on the difference of 
■02 and ipi. Thick grey lines are used in the following Figures to show the 
boundaries of Oj. 

Subsection 4.5. Avoiding the observer. Here, we find optimal paths 
for a robot navigating in a region containing stationary obstacles and a 
stationary enemy observer. The robot strives to be minimally exposed to 
the enemy at the same time making sure to avoid the obstacles and stay 
within the specified fuel budget. Observability by the enemy becomes the 
primary cost {Jo) to be minimized. As long as the secondary running cost K\ 
remains positive, our numerical method can solve the PDE (|17p efficiently 
even with with Kq = observability cost in f2j. However, this leads to 
infinitely many possible optimal paths since any portion of the path inside 
Cli would not contribute anything to Jq. To remove this arbitrariness and 
non- uniqueness, we assigned a small non-zero observability cost Kq = 0.1 on 
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Table 2 shows the other numerical parameters used in this experiment. 



primary: Observability, secondary: Path length 




[; , , , 
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Figure 1 1 . Optimal path for a robot navigating to minimize 
the exposure to a static enemy observer. 



Grid points in system-domain 


251 x 251 


Grid points along secondary budget 


301 


Speed of the vehicle 


1 


Primary cost (Observability cost) 


10 in the visible region 




0.1 in the non- visible region 


Secondary cost (Fuel rate) 


1 


Starting point for optimal path 


(0.1, 0.1) 


Target point for optimal path 


(1.0, 1.0) 


Observer location 


(0.15, 0) 



Table 2. Numerical parameters for the exposure- 
minimization example (Figure [TT]) . 



Figure [TT] shows the non- visible region with its boundary as thick solid 
lines. The small rectangles represent obstacles. Optimal paths correspond- 
ing to different bounds B\ on fuel budget are also plotted. Since the speed of 
motion is constant, the running costs are piecewise constant, and the obsta- 
cles are polygonal, it is easy to prove that all constraint-optimal paths are 
piecewise linear. As expected, the primary-optimal paths creep along the 
boundary of non-visible region. In this example, the target lies in Qf, the 
robot therefore follows the straight line path after it enters that component 

of f2 j. 
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Subsection 4.6. Minimizing fuel consumption/path length, con- 
strained by enemy observability. Here, we find the optimal path for 
a robot with the same two objectives as in the previous section. But now 
the path length (or fuel consumption) is the primary cost and the enemy ob- 
servability is secondary. The numerical parameters are shown in Table (3). 



primary: Fuel budget, secondary: Observability 




.25 .50 .75 1.0 



Figure 12. Optimal path for a robot minimizing its fuel 
budget/path length under the maximum enemy exposure 
constraint. The enemy exposure budget of 4.62 is large 
enough to allow following the primary-optimal trajectory 
(thick line). 



Grid points in system-domain 


301 x 301 


Grid points along secondary budget 


501 


Speed of the vehicle 


1 


Primary cost (Fuel rate/Path length) 


1 


Secondary cost (Observability) 


1 in the non-visible region 




5 in the visible region 


Starting point for optimal path 


(0.1, 0.1) 


Target point for optimal path 


(1.0, 1.0) 


Observer location 


(0.85, 0) 



Table 3. Numerical parameters for path length minimiza- 
tion constrained by exposure (Figure [12]) . 



Figure [T2l shows the optimal paths for this example. The secondary-optimal 
path shown as a dashed line is the least exposed to the enemy. As expected, 
the primary-optimal paths become shorter as we relax the constraint (or, 
alternatively, increase the "budget") for the maximum observability. 
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Subsection 4.7. Striving to be observed. Given a stationary friendly 
observer, the goal in many applications is to minimize the total time outside 
of direct visibility while moving to the target. This is more or less the op- 
posite of the problem considered in section 14.51 The total fuel available (or 
the maximum path length) is still treated as a constraint. The numerical 
parameters are shown in table (4). 

Figure [13] plots the optimal paths corresponding to different bounds on 
path length. When the bound on path length is tight, the robot has no 
option but to navigate through Oj. As we relax this bound, the robot finds 
a path which is always exposed to the observer. 

primary: non-observability, secondary: Path length 

1.00 



.75 
.50 



.25 


.25 .50 .75 1.00 



Figure 13. Optimal path for a robot in the presence of a 
friendly observer. 



Grid points in system-domain 


201 x 201 


Grid points along secondary budget 


301 


Speed of the vehicle 


1 


Primary cost (Non-observability) 


5 in non-visible region 




1 in visible region 


Secondary cost (Fuel rate) 


1 


Starting point for optimal path 


(0.1, 0.1) 


Target point for optimal path 


(1.0, 1.0) 


Observer location 


(0.85, 0) 



Table 4. Numerical parameters for minimizing the non- 



exposure constrained by path length (Figure [13]) . 
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Subsection 4.8. Path length minimization subject to two integral 
constraints. In this last example, we consider a problem of finding constrained- 
optimal paths in the presence of obstacles and two observers. The goal 
is to minimize the path-length subject to constraints on the amounts of 
time the robot can be visible to the enemy observer and invisible to the 
friendly observer. Given two secondary costs, the numerical domain is four- 
dimensional. As explained in section 13.51 we first solve the PDE (j24[) on 
Cl x [0, B<2\ to find the feasibility surface. We then march in the direction b\ 
to solve for the value function w. Two constrained-optimal trajectories are 
shown in Figure HH 



Grid points in system-domain 


101 x 101 


Grid points along each secondary budget 


301 


Speed of the vehicle 


1 


Primary cost 


path length 


Secondary cost 1 


visibility to enemy 


Secondary cost 2 


invisibility to friend 



Visibility cost values: 1 and 10 for both enemy and observer. 

Table 5. Numerical parameters for the two-secondary-costs 
example (Figure [T4"|) . 



Subsection 4.9. Discussion of computational complexity. Since our 
semi-Lagrangian discretization of the PDE (|17p is explicitly causal, the com- 
putational complexity of the methods is O(M), where M is the number of 
gridpoint on Q. e . 

A careful restriction of the problem to a feasible subdomain yields an 
efficient numerical method. For the case r = 1, computations on a 3- 
dimensional grid are quite fast on an average laptop. We have used Dell 
Inspiron 1505 laptop with 2 GHz Intel Centrino processor, and 1 GB RAM. 
On a grid with M = 201 3 gridpoints our instrumented (and unoptimized) 
code took 21, 43, and 39 seconds for Algorithms 1, 2, and 3 of section [3] 
respectively. For r = 2 and M = 101 4 , Algorithm 2 runs in less than 30 
minutes on the same laptop. 

In principle, only two &i-slices of the grid are needed in RAM to enable 
efficient marching. However, our current implementation allocates the entire 
grid. As a result, the last example (involving r = 2 and M = 101 2 x 301 2 
gridpoints and B\ = 11, Bi = 9) was computed on a machine with 64 
GB RAM though the memory footprint of the program is ~ 10 GB. The 
initialization (by solving for w\ on a 101 2 x 301 grid) took 17 seconds. The 
Algorithms 1, 2, and 3 then took 27, 67, and 55 minutes respectively. 

We note that the computational time is heavily dependent on the values 
of B\ and B2. There are two reasons for this. First, the tighter constraints 
make a larger part of Ct e non-feasible, resulting in a big reduction in com- 
putational cost for all three algorithms. Second, in Algorithms 2 and 3, r a 
is usually dependent on ratios between h and A6j's. This also influences 
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the number of (n + r)-dimensional cells traversed before reaching a suitable 
interpolation point (x,b); see section [3761 E.g., for r = 1, when Abi « h 
the interpolation is performed after traversing a single cell. Decreasing Bi 
while holding constant the number of gridpoints in that direction decreases 
Abi proportionally. To illustrate this point, the problem of section 14.81 on 
the same grid (M = 101 2 x 301 2 ) , but with B 1 = 4 and B 2 = 9 is much 
less computationally expensive: Algorithms 1,2, and 3 now take only 20, 27, 
and 26 minutes respectively on the same computer. 

In comparing these execution times to those reported in [25J, it is impor- 
tant to keep in mind that Mitchell and Sastry have used an upwind finite- 
difference discretization of the Eikonal PDE (i.e., isotropic control problems 
only). In contrast, our semi-Lagrangian implementation is also suitable for 
much more general (anisotropic and/or non-small-time-controllable) prob- 
lems, including those, where the minimization in (|17p cannot be performed 
analytically. 

We also note that our current implementation is simple but non-optimal 
since the grid is stored as a multidimensional array and all the domain- 
reduction procedures described in section 13.51 are performed only after the 
memory for the grid has been allocated. The "excluded" gridpoints are 
marked (to avoid computing W), but still take some computational time 
(in enumerations, input/output operations, etc.) A more efficient imple- 
mentation could be clearly built to allocate the memory for non-excluded 
gridpoints only, but this would require using non-array data structures to 
represent the grid. 

Section 5. Conclusions. 

We have introduced a new numerical method for multiobjective optimal 
control and single-objective optimal control in the presence of integral con- 
straints. Our approach is based on expanding the state space to include 
constraint-budgets and then solving an augmented PDE, whose explicit 
causality allows for a non-iterative (marching) numerical method. We have 
also shown the connection between the integral-constrained single-objective 
problem and the task of finding all Pareto-optimal controls. Our method 
was illustrated with a number of test-problems for two-dimensional optimal 
control with one and two additional integral constraints. We have used a 
flight-path bad weather avoidance example introduced in [25] as well as sev- 
eral examples of optimal robotic navigation in the presence of friendly and 
adversarial observers. 

It is a commonly accepted practical rule that lower-dimensional compu- 
tations are much less expensive than the higher-dimensional ones (simply 
because the number of gridpoints grows exponentially with the dimension). 
However, this simple rule of thumb ignores more subtle issues: How many 
PDEs need to be solved on each domain? How many times do we need to 
solve each PDE? Does the lower-dimensional approach adequately capture 
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the high-dimensional picture? Which PDE can be solved with a more effi- 
cient numerical method? 

These questions reflect the difference between our approach and the prior 
method by Mitchell and Sastry [25]. Their method is based on solving a 
system of (r + 2) PDEs (all but one of them linear, the remaining non-linear 
PDE with monotone causality, enabling a Dijkstra-like numerical method) 
on O C R n but with an r-dimensional parameter space, which requires 
solving this system repeatedly 0(2 r ) times. Our approach leads to a single 
PDE with explicit causality (enabling time-like marching) but on a (n + r)- 
dimensional domain. Even more importantly, our approach allows recover- 
ing the entire relevant part of Pareto front, including the non-convex parts 
of it, which are inaccessible using the weighted sums method employed in 
[25]. 

The explicit causality of the augmented PDE and a careful restriction of 
the problem to a feasible subdomain yield an efficient numerical method. 
However, the memory requirements of our method are more extensive since 
at least two &i-slices of the grid have to be kept in memory at all times for 
efficient marching. This is an (n + r — l)-dimensional grid, in contrast with 
an n-dimensional grid used by the method in [25]. Another disadvantage 
of our approach is the fact that the local truncation errors are 0(h) rather 
than 0(h). As a result, the quality of reconstruction of optimal controls and 
trajectories also depends on A&j's, whereas the method in [25] can provide 
a good trajectory reconstruction for each A regardless of coarseness of the 
mesh imposed on A. 

In the future we would like to build semi-Lagrangian and Eulerian higher- 
order accurate methods based on our approach. We also intend to explore 
the use of adaptive grids and unstructured meshes (since the constrained- 
optimal controls can be sensitive to small changes in available budgets). It 
will also be relevant to experiment with the efficiency /accuracy implications 
of the choice of marching direction - any other b% can be chosen since all 
Ki's are assumed to be positive, but the time r a in the semi-Lagrangian 
discretization is currently based on K\. 

Several other natural extensions should be possible without breaking the 
explicit causality of the augmented PDE. We intend to extend our method 
to 

(1) problems, where Kfs don't have to be positive for i > 1; 

(2) problems with running costs (and exit costs) dependent on the bud- 
gets still remaining; 

(3) optimal stochastic control subject to integral constraints; 

(4) differential games subject to integral constraints. 

The theoretical framework for considering secondary costs of varying sign 
has been developed in [37]. We believe that as long as Ki > holds at least 
for one i > 1, the computational efficiency will not be adversely affected. If 
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all secondary ZQ's (but not Kq) are allowed to change sign, this will require 
a method based on a more subtle monotone causality in w. 
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primary: Shortest path, b1 : Observability by Enemy, b2: Non-observability by Friend 
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Figure 14. A two-secondary costs example: fuel-optimality 

under constraints on visibility by enemy and non- visibility by 
friend. Dotted lines show the boundaries of visibility for both 
observers. Primary and secondary-optimal trajectories (top) 
and constrained-optimal trajectories (bottom). 



